Objekterkennung und Tracking á la CbTRAM

Der Ansatz aus Zinner et al. (2008) könnte interessant sein, um die Objektdefinition zu verbessern. Der Ansatz basiert auf drei Stadien. Interessant sind insbesondere die Definitionen der ersten beiden Stadien und die Verwendung der Bewegungsfelder.

Doch zunächst benötigen wir einige Pakete.

In [1]:
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

import datetime as dt

from scipy import ndimage as ndi

import sys
sys.path.append("/vols/talos/home/stephan/utils/tracking/")
import tracking_common as tc

from analysis_tools import optical_flow as oflow
import l15_msevi.msevi as msv

sys.path.append("/vols/talos/home/stephan/test/msevi_std2hires/")
import msevi_std2hires as msh
/vols/talos/local/anaconda2-5.0.0/lib/python2.7/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
/vols/talos/local/anaconda2-5.0.0/lib/python2.7/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools

Stadium 1: konvektive Auslösung

Der Ansatz für das erste Stadium geht vom HRV-Kanal aus. Die Grundidee ist, Wolken zu finden, die wachsen. Dazu wird das HRV-Bild eines Zeitpunktes t mit dem des Zeitpunktes t-1 verglichen. Das Bild des Zeitpunktes t-1 wird mittels eines Flussfeldes verschoben und dann die Differenz ΔHRV = HRV$_\mathrm{t}$ – HRV$_\mathrm{t-1}$ gebildet. Das Flussfeld wird dabei NICHT von den HRV-Bildern abgeleitet sondern aus den gelätteten IR10.8-Bildern der Zeisschritte t-2 und t-1, da die Bewegung im Flussfeld aus den HRV-Bildern nicht nur aus der Advektion sondern auch aus den Entwicklungen herrührt, die analysiert werden sollen.

Das ΔHRV wird mit einem ΔIR-Feld verglichen, dass so ähnlich wie das ΔHRV-Feld abgeleitet wird, nur dass das Bewegungsfeld dafür aus dem WV062-Kanal stammt. Da sich nicht alle Wolken im HRV konvektiv sind und sich auch nicht alle Cumuli humiles zu Cumolonimbi entwicklen, werden nur Wolken betrachtete, deren Flächenzunahme durch eine eine Abkühlung im IR-Feld gestützt wird. Dafür wird ein Qualitätsindikator aus den beiden Feldern ΔHRV und ΔIR gebildete. Dafür weden sie zunächst dann normalisiert und dann multipliziert um das Qualitätsfeld -ΔHRV $\cdot$ ΔIR zu erhalten, das Bereiche mit Wolkenwachstum anzeigt.

Um das auszuprobieren laden wir zunächst die Satellitendaten für ein Beispiel. Da das Verfahren vom HRV-Kanal ausgeht, müssen die Auflösungen angepasst werden. Damit keine Informationen verloren gehen, erhöhen wir die Auflösung der Kanäle mit Standardauflösung. Das sollte mit Hartwigs Werkzeug gehen.

In [68]:
t = t0
s = msv.MSevi(time=t,chan_list=sat_fields)
s.rad2bt()
s.rad2refl()
Region suggests use of hdf file
In [155]:
t2 = dt.datetime(2012,7,28,11,50) # t-2
t1 = dt.datetime(2012,7,28,11,55) # t-1
t0 = dt.datetime(2012,7,28,12,0) # t0

sat_fields = ['HRV','IR_108','WV_062']

region = ((216, 456), (1676, 2076))

sat_data = {f:[] for f in sat_fields}

for t in [t2,t1,t0]:
    s = msv.MSevi(time=t,chan_list=sat_fields,region=region)
    s.rad2bt()
    s.rad2refl()

    
    # Auflösung erhöhen
    for field in sat_fields[1:]:
               
        im = msh.upsample_image(s.bt[field],1,0)
        sat_data[field].append(im)
    
    sat_data['HRV'].append(s.ref['HRV'])
Region suggests use of hdf file
Region suggests use of hdf file
Region suggests use of hdf file

Als nächstes berechnen wir die benötigen Bewegungsfelder. Im CbTRAM-Artrikel wird das IR-Feld mit einem gleitenden Mittelwertfilter über 25 km geglättet, das entspricht also etwa 7 Standardpixeln oder 25 HRV-Pixeln.

In [156]:
ir_smooth = [ndi.filters.uniform_filter(ir_field,25) for ir_field in sat_data['IR_108']]
wv_smooth = [ndi.filters.uniform_filter(wv_field,25) for wv_field in sat_data['WV_062']]
In [157]:
fig,ax = plt.subplots(2,1,figsize = (16,20))
ax[0].imshow(sat_data['IR_108'][0],vmin=210,vmax=300,cmap='gray_r')
ax[0].set_title("Original")
ax[1].imshow(ir_smooth[0],vmin=210,vmax=300,cmap='gray_r')
ax[1].set_title(u"geglättet")
Out[157]:
Text(0.5,1,u'gegl\xe4ttet')

Nun können wir das Bewegungsfeld ermitteln.

In [158]:
ir_movement = tc.oft.calculate_optical_flow([sat_data['IR_108'][0],sat_data['IR_108'][1]],'farnebaeck')

Mit diesem Flussfeld, wird jetzt das HRV-Bild des Zeitpunktes t-1 auf t verschoben.

In [159]:
hrv_moved = oflow.morph_trans_opt_flow(sat_data['HRV'][-2],ir_movement[0])

Als nächstes berechnen wir die Differenz zwischen tatsächlichem und geschätzten Bild, um Bereiche mit Wolkenwachstum zu ermitteln.

In [161]:
DHRV = sat_data['HRV'][-1] - hrv_moved
In [162]:
def colourbar(mappable):
    ax = mappable.axes
    fig = ax.figure
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    return fig.colorbar(mappable, cax=cax)
In [163]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
d_plot = ax.imshow(DHRV,vmin=-0.5,vmax=0.5,cmap='RdBu_r')
ax.set_title(r"$\Delta$HRV")
colourbar(d_plot)
Out[163]:
<matplotlib.colorbar.Colorbar at 0x7fa996b75910>

Alle Bereiche mit einem positiven Wert weisen auf Wolkenwachstum hin. Um sicher zu gehen, betrachten wir die IR-Felder.

In [164]:
#wv_hires = [upsample_image(wv) for wv in sat_data['WV_062']]
wv_smooth = [ndi.filters.uniform_filter(wv_field,25) for wv_field in sat_data['WV_062']]
wv_movement = tc.oft.calculate_optical_flow([sat_data['WV_062'][0],sat_data['WV_062'][1]],'farnebaeck')
ir_moved = oflow.morph_trans_opt_flow(sat_data['IR_108'][-2],wv_movement[0])
In [165]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
d_plot = ax.imshow(ir_moved,vmin=210,vmax=300,cmap='gray_r')
colourbar(d_plot)
Out[165]:
<matplotlib.colorbar.Colorbar at 0x7fa994405f90>
In [166]:
DIR = sat_data['IR_108'][-1] - ir_moved
In [168]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
d_plot = ax.imshow(DIR,vmin=-10,vmax=10,cmap='RdBu_r')
ax.set_title(r"$\Delta$IR")
colourbar(d_plot)
Out[168]:
<matplotlib.colorbar.Colorbar at 0x7fa9940e8950>

Hier deuten alle Bereiche mit einer Abkühlung auf Wolkenwachstum hin.

Als nächstes werden beide Felder normalisiert.

In [300]:
DIR_norm = tc.scale_array_min_max(DIR,new_max=1,new_min=-1)
DHRV_norm = tc.scale_array_min_max(DHRV,new_max=1,new_min=-1)
In [303]:
fig,ax = plt.subplots(2,1,figsize=(16,20))
hrv_plot = ax[0].imshow(DHRV,vmin=-1,vmax=1,cmap='RdBu_r')
colourbar(hrv_plot)
ax[0].set_title(r"$\Delta$HRV")
hrv1_plot = ax[1].imshow(DHRV_norm,vmin=-1,vmax=1,cmap='RdBu_r')
colourbar(hrv1_plot)
ax[1].set_title(r"$\Delta$HRV, normalisiert")
Out[303]:
Text(0.5,1,'$\\Delta$HRV, normalisiert')

Daraus wird jetzt der Indikator für das Wolkenwachstum berechnet.

In [304]:
indicator = -DHRV_norm * DIR_norm
In [305]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
indd_plot = ax.imshow(indicator,vmin=-.5,vmax=.5,cmap='RdBu')
ax.set_title(r"-$\Delta$HRV $\cdot \Delta$IR")
colourbar(indd_plot)
Out[305]:
<matplotlib.colorbar.Colorbar at 0x7fa977336350>

Die blauen Bereiche sind Regionen mit Wolkenwachstum und die roten welche mit Wolkenauflösung.

Um wirklich interessante Objekte zu erhalten, muss ein Schwellwert auf das Feld angewendet werden. In Zinner et at. (2008) wird dreimal die Standardabweichung für die ΔHRV- und ΔIR-Felder verwendet. Wie der Schwellwert genau berechnet wird, wird nicht angegeben. Da die Felder im Bereich -1 bis 1 normalisiert werden, sollte, wenn sie normal verteilt sind, der Mittelwert 0 sein. Daher sollte $3\sigma$ oder $-3\sigma$ direkt als Schwellwert benutzt werden können. Für den HRV-Kanal kommt noch hinzu, dass alle Bereiche mit einer Reflektivität von wengier als 0,5 ausmaskiert werden.

In [330]:
threshold_DHRV = 3* np.std(DHRV_norm)
threshold_DIR = - 3* np.std(DIR_norm)
In [331]:
threshold_DHRV
Out[331]:
0.2984714807679782
In [332]:
threshold_DIR 
Out[332]:
-0.2362913088306438
In [333]:
threshold_indicator = -threshold_DHRV * threshold_DIR
threshold_indicator
Out[333]:
0.0705262168392859
In [343]:
DHRV_objects = np.ma.masked_less(DHRV_norm,threshold_DHRV)
DHRV_objects = np.ma.masked_less(DHRV_objects,0.5)
In [344]:
plt.imshow(DHRV_objects*1)
Out[344]:
<matplotlib.image.AxesImage at 0x7fa976e941d0>
In [345]:
DIR_objects = np.ma.masked_greater(DIR_norm,threshold_DIR)
plt.imshow(DIR_objects*1)
Out[345]:
<matplotlib.image.AxesImage at 0x7fa976e83b50>
In [357]:
hrv_mask = np.ma.masked_less(np.clip(sat_data['HRV'][-1],0,1),0.5)
level1_mask = np.ma.masked_less(indicator,threshold_indicator)

level1_mask = ~hrv_mask.mask & ~level1_mask.mask
In [359]:
fig,ax = plt.subplots(2,1,figsize=(16,20))
ax[0].imshow(~level1_mask,vmin=-1,vmax=1,cmap='gray')
ax[0].set_title(u"Maske mit Objekte im Stadium 1")
ax[1].imshow(np.clip(sat_data['HRV'][-1],0,1),vmin=0,vmax=1,cmap='gray')
ax[1].contour(~level1_mask*1,levels=1,colors='red')
ax[1].set_title(u"Objekte mit konvektiver Auslösung")
Out[359]:
Text(0.5,1,u'Objekte mit konvektiver Ausl\xf6sung')

Das sieht im Großen und Ganzen recht vernünftig aus, ist aber ziemlich auf die Randbereiche bereits existerender Konvektion fixiert.

Stadium 2: rasche Abkühlung

Dann versuchen wir mal Objekte des 2. Stadiums zu finden.

Dafür gehen wir genauso vor, wie zuvor, nur dass wir jetzt den WV062-Kanal benutzen. Der Fluss stammt auch von diesem Kanal, aber von den zwei Zeitschritten zuvor.

In [315]:
wv_moved = oflow.morph_trans_opt_flow(sat_data['WV_062'][1],wv_movement[0])

DWV = sat_data['WV_062'][-1] - wv_moved
In [316]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
wv_plot = ax.imshow(DWV,vmin=-5,vmax=5,cmap='binary')
colourbar(wv_plot)
ax.set_title(r"$\Delta$WV")
Out[316]:
Text(0.5,1,'$\\Delta$WV')

Die hellen Bereiche reflektieren Gebiete mit Abkühlung und die dunklen welche mit Erwärmung.

In [325]:
DWV_norm = tc.scale_array_min_max(DWV,new_max=1,new_min=-1)
#threshold_DWV = np.mean(DWV_norm) - 3 * np.std(DWV_norm)
threshold_DWV = - 3 * np.std(DWV_norm)

#level2_mask = np.ma.masked_greater(DWV_norm,-threshold_DWV)
level2_mask = np.ma.masked_greater(DWV_norm,threshold_DWV)
In [326]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
wv_plot = ax.imshow(DWV_norm,vmin=-1,vmax=1,cmap='binary')
colourbar(wv_plot)
ax.set_title(r"$\Delta$WV, normalisiert")
Out[326]:
Text(0.5,1,'$\\Delta$WV, normalisiert')
In [320]:
3 * np.std(DWV_norm)
Out[320]:
0.19226223906083953
In [327]:
plt.hist(DWV_norm.ravel(),bins=40)
Out[327]:
(array([3.00000e+00, 3.00000e+00, 1.00000e+00, 8.00000e+00, 1.30000e+01,
        1.40000e+01, 3.80000e+01, 7.00000e+01, 6.50000e+01, 1.33000e+02,
        1.30000e+02, 1.58000e+02, 2.13000e+02, 3.29000e+02, 4.94000e+02,
        7.30000e+02, 1.14000e+03, 2.15100e+03, 4.87000e+03, 1.07930e+04,
        2.71170e+04, 1.18960e+05, 4.24495e+05, 2.03333e+05, 4.36080e+04,
        1.51100e+04, 5.74800e+03, 2.20100e+03, 9.46000e+02, 4.96000e+02,
        2.55000e+02, 1.58000e+02, 8.00000e+01, 4.10000e+01, 3.30000e+01,
        2.20000e+01, 2.90000e+01, 7.00000e+00, 1.00000e+00, 4.00000e+00]),
 array([-1.  , -0.95, -0.9 , -0.85, -0.8 , -0.75, -0.7 , -0.65, -0.6 ,
        -0.55, -0.5 , -0.45, -0.4 , -0.35, -0.3 , -0.25, -0.2 , -0.15,
        -0.1 , -0.05,  0.  ,  0.05,  0.1 ,  0.15,  0.2 ,  0.25,  0.3 ,
         0.35,  0.4 ,  0.45,  0.5 ,  0.55,  0.6 ,  0.65,  0.7 ,  0.75,
         0.8 ,  0.85,  0.9 ,  0.95,  1.  ]),
 <a list of 40 Patch objects>)
In [361]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
wv_plot = ax.imshow(~level2_mask.mask*1,vmin=0,vmax=1,cmap='binary_r')
colourbar(wv_plot)
Out[361]:
<matplotlib.colorbar.Colorbar at 0x7fa975fe5450>
In [362]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
ax.imshow(np.clip(sat_data['HRV'][-1],0,1),vmin=0,vmax=1,cmap='gray')
ax.contour(~level1_mask*1,levels=1,colors='orange')
ax.contour(~level2_mask.mask*1,levels=1,colors='red')
Out[362]:
<matplotlib.contour.QuadContourSet at 0x7fa975f79750>
In [366]:
combined_masks = np.zeros_like(s.ref['HRV'])

combined_masks[np.where(level1_mask==True)] = 1
combined_masks[np.where(level2_mask.mask==False)] = 2
In [369]:
from matplotlib import colors

# make a color map of fixed colors
cmap = colors.ListedColormap([(0.9,0.9,0.9,0),'yellow','orange','red'])
bounds=[0,1,2,3,4]
norm = colors.BoundaryNorm(bounds, cmap.N)

fig,ax = plt.subplots(1,1,figsize=(16,10))
ax.imshow(np.clip(sat_data['HRV'][-1],0,1),vmin=0,vmax=1,cmap='gray')
level_plot = ax.imshow(combined_masks,cmap=cmap, norm=norm, alpha=0.4)
colourbar(level_plot)
Out[369]:
<matplotlib.colorbar.Colorbar at 0x7fa9763dcb90>

Das sieht im Großen ung Ganzen auch vernünftig aus.

In [368]:
import seaborn as sns
fig, ax = plt.subplots(1,1,figsize=(16,10))
levels = {1:u'konv. Auslösung', 2:'schn. Wachstum'}
for level in [1,2]:
    values = sat_data['IR_108'][-1][np.where(combined_masks==level)]
    
    # Draw the density plot
    sns.distplot(values, hist = False, kde = True,
                 kde_kws = {'linewidth': 3},
                 label = levels[level],ax=ax)
    
# Plot formatting
plt.legend(prop={'size': 16}, title = 'Stadium')
Out[368]:
<matplotlib.legend.Legend at 0x7fa975e4f550>